Experimentos com pêndulos

In [1]:
using Plots
using Images
using ImageDraw
using VideoIO
using Base64
using CSV
using DataFrames
using DifferentialEquations.OrdinaryDiffEq
using Optim
using StatsBase
In [2]:
include("scripts/DisplayMovie.jl")
include("scripts/ImageTracking.jl")
Out[2]:
Main.ImageTracking
In [3]:
filename = "img/pendulo_70cm_1_reduzido_curto.mp4"
vd = VideoIO.load(filename)
Out[3]:
(a vector displayed as a row to save space)
In [4]:
typeof(copy(vd[1]))
Out[4]:
Matrix{RGB{N0f8}} (alias for Array{RGB{Normed{UInt8, 8}}, 2})
In [5]:
println("Duração do vídeo: $(VideoIO.get_duration(filename)) s")
println("Número de quadros: $(length(vd))")
println("Número médio de quadros por segundo: $(length(vd)/VideoIO.get_duration(filename))")
Duração do vídeo: 5.086 s
Número de quadros: 150
Número médio de quadros por segundo: 29.49272512780181
In [6]:
# video completo
filename_completo = "img/pendulo_70cm_1_reduzido.mov"
vd_completo = VideoIO.load(filename_completo)
println("""Duração do vídeo: $(VideoIO.get_duration(filename_completo)) s""")
println("Número de quadros: $(length(vd_completo))")
println("""Número médio de quadros por segundo: $(length(vd_completo)/VideoIO.get_duration(filename_completo))""")
Duração do vídeo: 16.45 s
Número de quadros: 493
Número médio de quadros por segundo: 29.969604863221885
In [7]:
vd[1]
Out[7]:
In [8]:
DisplayMovie.display_movie("img/pendulo_70cm_1_reduzido_curto.mp4")
In [9]:
background = convert.(RGB{Float16}, vd[1])
Out[9]:
In [10]:
convert.(RGB{Float16}, vd[54])
Out[10]:
In [11]:
convert.(Gray{Float16}, abs.(convert.(RGB{Float16}, vd[54]) - background))
Out[11]:
In [12]:
threshold = 0.1
mask = convert.(
    Gray{Float16},
    abs.(
        convert.(RGB{Float16}, vd[54]) - 
        background
    )
) .> threshold
Out[12]:
480×854 BitMatrix:
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  1  1  1  1  1  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  1  1  1  1  1  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     1  0  1  1  1  1  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     1  0  1  1  1  1  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  1  0  0  0  0  1  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     1  1  0  0  0  0  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     1  0  1  0  0  0  0  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     1  0  1  0  0  0  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     1  1  0  1  0  0  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  1  1  0  0  1  1  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  1  0  0  1  1  1  1  1  1  1  1
 0  0  0  0  0  0  0  0  0  0  0  0  0     1  1  0  0  1  0  0  0  1  1  1  1
 ⋮              ⋮              ⋮        ⋱           ⋮              ⋮        
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  0  0  0     0  0  0  0  0  0  0  0  0  0  0  0
In [13]:
convert.(Gray{Float16}, mask)
Out[13]:
In [14]:
labeled_components = label_components(mask)
Out[14]:
480×854 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0     249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0     249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0     249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0     249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0  249  249  249  249
 ⋮              ⋮              ⋮        ⋱                   ⋮            
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0  …    0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0  …    0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
In [15]:
@which label_components(mask)
Out[15]:
In [16]:
?label_components
search: label_components label_components! labeled_components

Out[16]:
label = label_components(tf, [connectivity])
label = label_components(tf, [region])

Find the connected components in a binary array tf. There are two forms that connectivity can take:

  • It can be a boolean array of the same dimensionality as tf, of size 1 or 3 along each dimension. Each entry in the array determines whether a given neighbor is used for connectivity analyses. For example, connectivity = trues(3,3) would use 8-connectivity and test all pixels that touch the current one, even the corners.
  • You can provide a list indicating which dimensions are used to determine connectivity. For example, region = [1,3] would not test neighbors along dimension 2 for connectivity. This corresponds to just the nearest neighbors, i.e., 4-connectivity in 2d and 6-connectivity in 3d.

The default is region = 1:ndims(A).

The output label is an integer array, where 0 is used for background pixels, and each connected region gets a different integer index.

In [17]:
colorset = distinguishable_colors(maximum(labeled_components)+1, colorant"black", dropseed = false)
Out[17]:
In [18]:
labeled_components
Out[18]:
480×854 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0     249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0     249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0     249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0  …  249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0     249  249  249  249  249  249  249
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0  249  249  249  249
 ⋮              ⋮              ⋮        ⋱                   ⋮            
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0  …    0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0  …    0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
 0  0  0  0  0  0  0  0  0  0  0  0  0       0    0    0    0    0    0    0
In [19]:
coloredmask = map( n -> colorset[n+1], labeled_components)
Out[19]:
In [20]:
maximum(labeled_components)
Out[20]:
261
In [21]:
components_location = ImageTracking.locate_components(labeled_components)
min_area = 1
filter!(p -> p.area ≥ min_area, components_location)
length(components_location)
Out[21]:
261
In [22]:
components_location = ImageTracking.locate_components(labeled_components)
min_area = 10
filter!(p -> p.area ≥ min_area, components_location)
length(components_location)
Out[22]:
18
In [23]:
vd54 = copy(vd[54])
for p in components_location
    draw!(vd54, Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
end
vd54
Out[23]:
In [24]:
components_location = ImageTracking.locate_components(labeled_components)
min_area = 100
filter!(p -> p.area ≥ min_area, components_location)
Out[24]:
3-element Vector{Main.ImageTracking.Blob}:
 Blob Object
  Centroid coordinates: (x, y) = (670.1450079239303, 387.1695721077655)
  Span values: (xspan, yspan) = (640:701, 360:413)
  Occupied: area = 2524

 Blob Object
  Centroid coordinates: (x, y) = (842.4875, 20.16875)
  Span values: (xspan, yspan) = (837:849, 1:48)
  Occupied: area = 160

 Blob Object
  Centroid coordinates: (x, y) = (850.0762711864406, 7.669491525423729)
  Span values: (xspan, yspan) = (844:854, 1:17)
  Occupied: area = 118
In [25]:
vd54 = copy(vd[54])
for p in components_location
    draw!(vd54, Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
end
vd54
Out[25]:
In [26]:
components_location = ImageTracking.locate_components(labeled_components)
min_area = 200
filter!(p -> p.area ≥ min_area, components_location)
Out[26]:
1-element Vector{Main.ImageTracking.Blob}:
 Blob Object
  Centroid coordinates: (x, y) = (670.1450079239303, 387.1695721077655)
  Span values: (xspan, yspan) = (640:701, 360:413)
  Occupied: area = 2524
In [27]:
p = components_location[1]
draw(vd[54], Polygon(ImageTracking.RectanglePoints(p)), colorant"red")
Out[27]:
In [28]:
filename = "img/pendulo_70cm_1_reduzido_curto.mp4"
tracks_70cm1 = ImageTracking.find_tracks(filename)
Processing 150 frames ...100%|██████████████████████████| Time: 0:00:47
Out[28]:
1-element Vector{Main.ImageTracking.Track}:
 Tracked blob with framespan nspan = 13:150
In [29]:
ImageTracking.tracks_on_video(filename, tracks_70cm1)
┌ Info: Video saved with 150 frames and 1 tracks
â”” @ Main.ImageTracking /Users/rrosa/Documents/git_repositories/modelagem_matematica/notas_de_aula/scripts/ImageTracking.jl:150
In [30]:
DisplayMovie.display_movie("img/pendulo_70cm_1_reduzido_curto_tracked.mp4")
In [31]:
tracks_70cm1
Out[31]:
1-element Vector{Main.ImageTracking.Track}:
 Tracked blob with framespan nspan = 13:150
In [32]:
tracks_70cm1[1].nspan
Out[32]:
13:150
In [33]:
tracks_70cm1[1].path
Out[33]:
138-element Vector{Main.ImageTracking.Blob}:
 Blob Object
  Centroid coordinates: (x, y) = (848.7265500794913, 209.04451510333863)
  Span values: (xspan, yspan) = (839:854, 175:241)
  Occupied: area = 629

 Blob Object
  Centroid coordinates: (x, y) = (847.3719862227325, 212.89322617680827)
  Span values: (xspan, yspan) = (836:854, 176:248)
  Occupied: area = 871

 Blob Object
  Centroid coordinates: (x, y) = (846.1327913279133, 216.11020776874435)
  Span values: (xspan, yspan) = (833:854, 177:254)
  Occupied: area = 1107

 Blob Object
  Centroid coordinates: (x, y) = (844.8869047619048, 218.23511904761904)
  Span values: (xspan, yspan) = (830:854, 177:256)
  Occupied: area = 1344

 Blob Object
  Centroid coordinates: (x, y) = (843.5773130544993, 220.33523447401774)
  Span values: (xspan, yspan) = (828:854, 177:258)
  Occupied: area = 1578

 Blob Object
  Centroid coordinates: (x, y) = (842.7289227166276, 221.2980093676815)
  Span values: (xspan, yspan) = (826:854, 177:258)
  Occupied: area = 1708

 Blob Object
  Centroid coordinates: (x, y) = (842.6964285714286, 220.60059523809525)
  Span values: (xspan, yspan) = (826:854, 178:258)
  Occupied: area = 1680

 Blob Object
  Centroid coordinates: (x, y) = (842.8163514338011, 220.41122635753507)
  Span values: (xspan, yspan) = (826:854, 178:258)
  Occupied: area = 1639

 Blob Object
  Centroid coordinates: (x, y) = (842.8624923265808, 220.86126457949663)
  Span values: (xspan, yspan) = (826:854, 179:259)
  Occupied: area = 1629

 Blob Object
  Centroid coordinates: (x, y) = (842.980625, 221.265625)
  Span values: (xspan, yspan) = (826:854, 179:258)
  Occupied: area = 1600

 Blob Object
  Centroid coordinates: (x, y) = (842.8601615910503, 221.1224362958359)
  Span values: (xspan, yspan) = (826:854, 179:258)
  Occupied: area = 1609

 Blob Object
  Centroid coordinates: (x, y) = (842.5529622980251, 221.25852782764812)
  Span values: (xspan, yspan) = (825:854, 179:259)
  Occupied: area = 1671

 Blob Object
  Centroid coordinates: (x, y) = (842.5005966587112, 221.49224343675417)
  Span values: (xspan, yspan) = (825:854, 179:259)
  Occupied: area = 1676

 â‹®
 Blob Object
  Centroid coordinates: (x, y) = (537.5756442227764, 423.214463840399)
  Span values: (xspan, yspan) = (507:569, 399:446)
  Occupied: area = 2406

 Blob Object
  Centroid coordinates: (x, y) = (586.1208315288003, 414.7986141186661)
  Span values: (xspan, yspan) = (556:617, 391:438)
  Occupied: area = 2309

 Blob Object
  Centroid coordinates: (x, y) = (631.0521201413428, 402.50132508833923)
  Span values: (xspan, yspan) = (602:661, 379:426)
  Occupied: area = 2264

 Blob Object
  Centroid coordinates: (x, y) = (671.4093648585999, 387.42280945757994)
  Span values: (xspan, yspan) = (644:701, 364:411)
  Occupied: area = 2157

 Blob Object
  Centroid coordinates: (x, y) = (707.1664247701983, 370.4446057087566)
  Span values: (xspan, yspan) = (681:735, 348:394)
  Occupied: area = 2067

 Blob Object
  Centroid coordinates: (x, y) = (737.0905852417303, 352.37150127226465)
  Span values: (xspan, yspan) = (712:764, 330:376)
  Occupied: area = 1965

 Blob Object
  Centroid coordinates: (x, y) = (762.6097691894794, 334.12238325281805)
  Span values: (xspan, yspan) = (738:788, 312:358)
  Occupied: area = 1863

 Blob Object
  Centroid coordinates: (x, y) = (783.1237172177879, 317.3369441277081)
  Span values: (xspan, yspan) = (759:808, 296:340)
  Occupied: area = 1754

 Blob Object
  Centroid coordinates: (x, y) = (798.8568935427575, 302.46364165212333)
  Span values: (xspan, yspan) = (775:823, 281:325)
  Occupied: area = 1719

 Blob Object
  Centroid coordinates: (x, y) = (810.6108343711084, 290.052303860523)
  Span values: (xspan, yspan) = (787:834, 269:311)
  Occupied: area = 1606

 Blob Object
  Centroid coordinates: (x, y) = (818.3419857235561, 280.7988319273199)
  Span values: (xspan, yspan) = (796:841, 260:302)
  Occupied: area = 1541

 Blob Object
  Centroid coordinates: (x, y) = (822.5608930987821, 274.98240866035184)
  Span values: (xspan, yspan) = (800:845, 255:295)
  Occupied: area = 1478
In [34]:
tracks_70cm1[1].path[1]
Out[34]:
Blob Object
  Centroid coordinates: (x, y) = (848.7265500794913, 209.04451510333863)
  Span values: (xspan, yspan) = (839:854, 175:241)
  Occupied: area = 629
In [35]:
plt = plot(title = "Paths of $(length(tracks_70cm1)) track(s)", titlefont = 10,
    xlims = (1,854), ylims = (1, 480),
    size = (854, 480), aspect = :equal, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
    scatter!(plt, getfield.(tr.path, :x), 481 .- getfield.(tr.path, :y), label="track $n")
end
display(plt)
In [36]:
plt = plot(title = "Coordinate x of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topleft)
for (n, tr) in enumerate(tracks_70cm1)
    scatter!(plt, tr.nspan, getfield.(tr.path, :x), label="track $n")
end
display(plt)
In [37]:
plt = plot(title = "Coordinate y of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
    scatter!(plt, tr.nspan, - getfield.(tr.path, :y), label="track $n")
end
display(plt)
In [38]:
size(vd[1])
Out[38]:
(480, 854)
In [39]:
data_x = getfield.(tracks_70cm1[1].path, :x)
data_y = size(vd[1], 1) + 1 .- getfield.(tracks_70cm1[1].path, :y)
nothing
In [40]:
n_minima = Int[]
for j in 40:length(data_x)-1
    if data_x[j-1] > data_x[j] < data_x[j+1]
        push!(n_minima, j + first(tracks_70cm1[1].nspan) - 1)
    end
end
n_minima
Out[40]:
2-element Vector{Int64}:
  72
 125
In [41]:
plt = plot(title = "Coordinate x of $(length(tracks_70cm1)) tracks", titlefont = 10, legend=:topleft)
for (n, tr) in enumerate(tracks_70cm1)
    scatter!(plt, tr.nspan, getfield.(tr.path, :x), label="track $n")
end
vline!(plt, n_minima, label = "minima")
display(plt)
In [42]:
(n_minima[2:end] - n_minima[1:end-1]) / 29.5
Out[42]:
1-element Vector{Float64}:
 1.7966101694915255
In [43]:
g = 9.8 # m/s
l = 0.7 # m
T_p = 2π * √(l / g)
Out[43]:
1.679251908362714
In [44]:
function medias(a::Real, b::Real)
    (a ≥ 0 && b ≥ 0) || throw(ArgumentError("arguments must be nonnegative"))
    ma = (a + b)/2
    mg = sqrt(a * b)
    return ma, mg
end

function agm(a::Real, b::Real; tol::Real = 1e-10, maxiter::Int = 100)
    tol > 0 || throw(ArgumentError("tolerance must be positive"))
    maxiter > 0 || throw(ArgumentError("maximum number of iterations must be positive"))
    y1, y2 = a, b
    n = 0
    while abs(y1 - y2) > tol && n < maxiter
        y1, y2 = medias(y1, y2)
        n += 1
    end
    return (y1 + y2)/2, abs(y1 - y2), n
end
Out[44]:
agm (generic function with 1 method)
In [45]:
θ₀ = π / 4
T_p / agm(1, cos(θ₀ / 2))[1]
Out[45]:
1.7463772212095845

Escala

  • Vamos fazer ajustar um circunferência ao movimento do pêndulo, para descobrir a escala espacial.
In [108]:
error(u) = sum(
    abs2,
    (data_x[50:end] .- u[1]).^2 + (data_y[50:end] .- u[2]).^2 .- u[3]^2
    ) / length(data_x)
Out[108]:
error (generic function with 1 method)
In [109]:
center_point = Optim.optimize(error, [400.0, 800.0, 500])
Out[109]:
 * Status: success

 * Candidate solution
    Final objective value:     2.179085e+06

 * Found with
    Algorithm:     Nelder-Mead

 * Convergence measures
    √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    152
    f(x) calls:    280
In [110]:
Optim.minimum(center_point)
Out[110]:
2.1790845446313173e6
In [49]:
center_nx, center_ny, radius_n = Optim.minimizer(center_point)
Out[49]:
3-element Vector{Float64}:
 466.34634496986183
 558.3533645572189
 505.51516552719033
In [105]:
plt = plot(title = "Paths of $(length(tracks_70cm1)) track(s)", titlefont = 10,
    xlims = (1,854), ylims = (1, 600),
    size = (854, 600), aspect = :equal, legend=:topright)
for (n, tr) in enumerate(tracks_70cm1)
    scatter!(plt, data_x, data_y, label="track $n")
end
scatter!(plt, (center_nx, center_ny), markersize = 10)
plot!(plt, 1:854, nx -> center_ny - √(radius_n^2 - (nx - center_nx)^2), linestyle = :dash)
display(plt)
In [51]:
scale = radius_n / l # points per meter
Out[51]:
722.1645221817005
In [52]:
data_scaled_x = (data_x .- center_nx) / scale
data_scaled_y = (data_y .- center_ny) / scale
Out[52]:
138-element Vector{Float64}:
 -0.3965825942201275
 -0.4019120045625831
 -0.4063666426583695
 -0.40930906258292526
 -0.41221714704552687
 -0.41355032648607204
 -0.412584598998547
 -0.4123223749834595
 -0.4129455546165464
 -0.41350548300970785
 -0.41330720588619096
 -0.4134956553705288
 -0.413819287454254
  â‹®
 -0.6931492935783854
 -0.6814956475417342
 -0.6644672715240699
 -0.6435876586829319
 -0.6200774982868892
 -0.5950512004262687
 -0.569781171978575
 -0.5465379377714996
 -0.5259424889246751
 -0.5087561866204509
 -0.49594266331796594
 -0.48788851071379713
In [53]:
function dudt_pendulum!(du, u, p, t)
    θ, ω = u
    l, g, α = p
    du[1] = ω
    du[2] = - (g / l) * sin(θ) - α * ω
end
Out[53]:
dudt_pendulum! (generic function with 1 method)
In [54]:
g = 9.8 # 9,8 m/s^2
l = 0.7 # 70 cm
α = 0
p = [l, g, α]
θ₀ = π/6
ω₀ = 0.0
u₀ = [θ₀, ω₀]
tspan = (0.0, 5.0) # até 5 segundos
prob_pendulum = ODEProblem(dudt_pendulum!, uâ‚€, tspan, p)
Out[54]:
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 5.0)
u0: 2-element Vector{Float64}:
 0.5235987755982988
 0.0
In [55]:
sol = solve(prob_pendulum, Tsit5())
Out[55]:
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 32-element Vector{Float64}:
 0.0
 0.00014258482589892232
 0.0015684330848881453
 0.015826915674780374
 0.05520287938356151
 0.11869948603402829
 0.2002131852019463
 0.305914875953137
 0.42890164394056673
 0.5569187047083594
 0.7013506994917333
 0.872697562411765
 1.0408497431878225
 â‹®
 2.6895148788362473
 2.8890581891668203
 3.1312540989213207
 3.3360859425870872
 3.579811053395161
 3.7899434234522027
 4.033677808420313
 4.248777258478439
 4.489724412426027
 4.7129429751996135
 4.9487942885038185
 5.0
u: 32-element Vector{Vector{Float64}}:
 [0.5235987755982988, 0.0]
 [0.5235987044417862, -0.000998093740288535]
 [0.5235901656815024, -0.01097897701799843]
 [0.52272227807514, -0.1107323314966049]
 [0.512965861754178, -0.3840406673395072]
 [0.47498717025387077, -0.8072478892241957]
 [0.388974682500463, -1.2882714556007608]
 [0.22681168930963705, -1.74188417532362]
 [-0.0034489375342210073, -1.936781390502851]
 [-0.2416589452693402, -1.7139641695093224]
 [-0.44371540806352233, -1.0197244000033383]
 [-0.5224044417637903, 0.12912618721927982]
 [-0.40599595484510004, 1.214547948859705]
 â‹®
 [-0.468181398562196, 0.8593957023872523]
 [-0.19043787415408014, 1.8015157078868282]
 [0.26158347131616566, 1.6731851139158112]
 [0.5009792313597429, 0.5577249920462498]
 [0.4331262533935047, -1.0800627594148071]
 [0.10380928939438962, -1.8977491548843235]
 [-0.337721682376994, -1.4733201977988184]
 [-0.5219373906746808, -0.15523522675671586]
 [-0.3642262071170868, 1.384044759406596]
 [0.028820747076412908, 1.9341611787299822]
 [0.4182220223013869, 1.1572024377115584]
 [0.4697150849008621, 0.8484644252996695]
In [ ]:

In [56]:
display(plot(sol, label = ["θ(t)" "ω(t)"]))
display(plot(sol, vars = 1, label = "θ(t)"))
In [57]:
n0 = 34
plt = plot(sol, vars = 1, label = "θ(t)")
scatter!(plt, (0:length(data_x)-n0) ./ 29.5, data_scaled_x[n0:end], label="tracked")
display(plt)
In [58]:
df = CSV.read("data/pendulo_70cm_reduzido.csv", DataFrame)
Out[58]:

481 rows × 3 columns

nfnxny
Int64Float64Float64
113848.727209.073
214847.362212.843
315846.127216.006
416844.89218.281
517843.577220.369
618842.701221.178
719842.682220.634
820842.833220.468
921842.866220.873
1022843.008221.232
1123842.875221.118
1224842.552221.235
1325842.523221.488
1426842.712221.488
1527842.606221.49
1628842.406221.461
1729842.366221.539
1830842.363221.584
1931842.25221.574
2032842.034221.523
2133841.969221.37
2234841.981221.283
2335841.954221.323
2436842.037221.198
2537842.135221.06
2638842.189221.251
2739842.266221.351
2840842.162221.446
2941841.951221.716
3042841.715221.619
In [115]:
plt = plot(title = "Coordinate x", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, df.nx, label = nothing)
display(plt)

plt = plot(title = "Coordinate y", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, size(vd[1], 1) .+ 1 .- df.ny, label = nothing)
display(plt)
In [114]:
data2_scaled_t = (df.nf[n0:end] .- df.nf[n0]) / 29.97
data2_scaled_x = (df.nx[n0:end] .- center_nx) / scale
data2_scaled_y = (size(vd[1], 1) .+ 1 .- df.ny[n0:end] .- center_ny ) / scale
display(scatter(data2_scaled_t, data2_scaled_x))
display(scatter(data2_scaled_t, data2_scaled_y))
In [62]:
data_x2 = getfield.(tracks2[1].path, :x)
n_minima2 = Int[]
for j in 40:length(data_x2)-1
    if data_x2[j-1] > data_x2[j] < data_x2[j+1]
        push!(n_minima2, j + first(tracks2[1].nspan) - 1)
    end
end
n_minima2
In [116]:
n_minima2 = Int[]
for j in 40:length(df.nx)-1
    if df.nx[j-1] > df.nx[j] < df.nx[j+1]
        push!(n_minima2, j + df.nf[1] - 1)
    end
end
n_minima2
Out[116]:
8-element Vector{Int64}:
  72
 125
 178
 231
 283
 336
 389
 441
In [117]:
(n_minima2[2:end] - n_minima2[1:end-1])/30
Out[117]:
7-element Vector{Float64}:
 1.7666666666666666
 1.7666666666666666
 1.7666666666666666
 1.7333333333333334
 1.7666666666666666
 1.7666666666666666
 1.7333333333333334
In [118]:
sum(n_minima2[2:end] - n_minima2[1:end-1])/7/30
Out[118]:
1.7571428571428571
In [122]:
plt = plot(title = "Coordinate x", titlefont = 10, legend=:topleft)
scatter!(plt, df.nf, df.nx, label = "data")
vline!(plt, n_minima2, label = "minima")
display(plt)
In [66]:
g = 9.8 # 9,8 m/s^2
l = 0.7 # 70 cm
α = 0
p = [l, g, α]
θ₀ = π/6
ω₀ = 0.0
u₀ = [θ₀, ω₀]
tspan2 = (0.0, 16.45)
prob_pendulum2 = ODEProblem(dudt_pendulum!, uâ‚€, tspan2, p)
Out[66]:
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 16.45)
u0: 2-element Vector{Float64}:
 0.5235987755982988
 0.0
In [67]:
VideoIO.get_duration("img/pendulo_70cm_1_reduzido.mov")
Out[67]:
16.45
In [68]:
θ₀ = π / 6
α = 0.0 #-0.1
remake(prob_pendulum2; u0 = [θ₀, 0.0], p = [l, g, α])
sol2 = solve(prob_pendulum2, Tsit5(); saveat = 1/29.87)
data2_scaled_x[1:10] .- sol2(data2_scaled_t)[1,1:10]
#length(sol[1,:])
#sol(data2_scaled_t)[1,:]
sol(data2_scaled_t)
Out[68]:
t: 448-element Vector{Float64}:
  0.0
  0.033366700033366704
  0.06673340006673341
  0.1001001001001001
  0.13346680013346682
  0.16683350016683351
  0.2002002002002002
  0.2335669002335669
  0.26693360026693363
  0.3003003003003003
  0.33366700033366703
  0.3670337003670337
  0.4004004004004004
  â‹®
 14.547881214547882
 14.58124791458125
 14.614614614614615
 14.647981314647982
 14.68134801468135
 14.714714714714715
 14.748081414748082
 14.78144811478145
 14.814814814814815
 14.848181514848182
 14.88154821488155
 14.914914914914915
u: 448-element Vector{Vector{Float64}}:
 [0.5235987755982988, 0.0]
 [0.5197064819391034, -0.2330414413613375]
 [0.5080822001687944, -0.4629302691579599]
 [0.488883698677749, -0.6865155884298596]
 [0.46237373952364935, -0.9006563907398047]
 [0.4289196227578753, -1.102237149997204]
 [0.38899141026001377, -1.288202512164673]
 [0.3431598607271923, -1.4556077279142619]
 [0.29209180679967006, -1.6016605695396418]
 [0.23654135225424, -1.7238352335082683]
 [0.17734352255250205, -1.819932863153441]
 [0.11540186183275535, -1.88814487747045]
 [0.05166798968064211, -1.9271578012356183]
 â‹®
 [26920.697182912332, 33556.85056390745]
 [27306.360007144147, 33989.425031612314]
 [27696.14382775103, 34426.118397211016]
 [28090.07784082385, 34866.95635621334]
 [28488.191345596126, 35311.96468510579]
 [28890.51374518916, 35761.16921555201]
 [29297.07454454987, 36214.59588839963]
 [29707.903351789308, 36672.27070221071]
 [30123.029879305846, 37134.21973321311]
 [30542.483937377587, 37600.4691636604]
 [30966.29544585134, 38071.04522386799]
 [31394.494424368553, 38545.97422357512]
In [69]:
plot(data2_scaled_t, sol(data2_scaled_t)[1,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")
Out[69]:
In [70]:
function ddu_pendulum!(ddu, du, u, p, t)
    θ = u[1]
    ω = du[1]
    l, g, α = p
    ddu[1] = - (g / l) * sin(θ) - α * ω
end
Out[70]:
ddu_pendulum! (generic function with 1 method)
In [71]:
θ₀ = π / 6
ω₀ = 0.0
α = 0.05
l = 0.74
p = [l, g, α]
prob_pendulum_2nd = SecondOrderODEProblem(ddu_pendulum!, [ω₀], [θ₀], tspan2, p)
sol_2nd = solve(prob_pendulum_2nd)
sol_2nd.retcode
Out[71]:
:Success
In [72]:
plot(data2_scaled_t, sol_2nd(data2_scaled_t)[2,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")
Out[72]:
In [137]:
prob_rmk = remake(prob_pendulum_2nd; u0 = [θ₀], p = [l, g, α])
sol = solve(prob_rmk)
type Array has no field x

Stacktrace:
  [1] getproperty(x::Vector{Float64}, f::Symbol)
    @ Base ./Base.jl:33
  [2] (::DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing})(du::Vector{Float64}, u::Vector{Float64}, p::Vector{Float64}, t::Float64)
    @ SciMLBase ~/.julia/packages/SciMLBase/UIp7W/src/scimlfunctions.jl:354
  [3] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{true}}, Tsit5, OrdinaryDiffEq.InterpolationData{DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/perform_step/low_order_rk_perform_step.jl:623
  [4] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{true}}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:default_set, :second_time), Tuple{Bool, Bool}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:455
  [5] #__solve#404
    @ ~/.julia/packages/OrdinaryDiffEq/2ZBfC/src/solve.jl:4 [inlined]
  [6] __solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{true}}, ::Nothing; default_set::Bool, kwargs::Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:second_time,), Tuple{Bool}}})
    @ DifferentialEquations ~/.julia/packages/DifferentialEquations/7WPQg/src/default_solve.jl:8
  [7] #__solve#57
    @ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:291 [inlined]
  [8] __solve
    @ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:278 [inlined]
  [9] #solve_call#42
    @ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:61 [inlined]
 [10] solve_call
    @ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:48 [inlined]
 [11] #solve_up#44
    @ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:90 [inlined]
 [12] solve_up
    @ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:78 [inlined]
 [13] #solve#43
    @ ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:73 [inlined]
 [14] solve(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, DynamicalODEFunction{true, ODEFunction{true, typeof(ddu_pendulum!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{true, SciMLBase.var"#225#227", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SecondOrderODEProblem{true}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/WKucm/src/solve.jl:68
 [15] top-level scope
    @ In[137]:2
 [16] eval
    @ ./boot.jl:360 [inlined]
 [17] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1116
In [94]:
function objective(β)
    g = 9.8
    ω₀ = 0.0
    θ₀, l, α = β
    p = [l, g, α]
    tspan2 = (0.0, 16.45)
    prob = remake(prob_pendulum_2nd, u0 = [θ₀, 0.0], p = [l, g, α])
    #sol = solve(prob)
    prob = SecondOrderODEProblem(ddu_pendulum!, [ω₀], [θ₀], tspan2, p)
    sol = solve(prob)
    return mse(data2_scaled_x, sol(data2_scaled_t)[2,:])
end
Out[94]:
objective (generic function with 1 method)
In [95]:
res = Optim.optimize(objective, [Ï€ / 3, 0.7, 0.0])
Out[95]:
 * Status: success

 * Candidate solution
    Final objective value:     1.313365e-03

 * Found with
    Algorithm:     Nelder-Mead

 * Convergence measures
    √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    69
    f(x) calls:    137
In [83]:
θ₀, l, α = Optim.minimizer(res)
Out[83]:
3-element Vector{Float64}:
 0.5557462725091609
 0.7421821761291346
 0.03752799438015029
In [76]:
π / θ₀
Out[76]:
5.652926180513441
In [77]:
θ₀ - π / 5.6529
Out[77]:
-2.57385107782504e-6
In [78]:
remake(prob_pendulum_2nd, u0 = [θ₀, 0.0], p = [l, g, α])
Out[78]:
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 16.45)
u0: 2-element Vector{Float64}:
 0.5557462725091609
 0.0
In [79]:
sol_2nd = solve(prob_pendulum_2nd)
sol_2nd.retcode
Out[79]:
:Success
In [80]:
plot(data2_scaled_t, sol_2nd(data2_scaled_t)[2,:], label = "ODE", ylims=(-1,1))
scatter!(data2_scaled_t, data2_scaled_x, label = "data")
Out[80]: